#!/usr/bin/perl
use strict;
# use warnings;

#####################################################

# my ($USAGE) =  "\n\n\t****** USAGE of preprocess.pl PROGRAM ******\n\n\n\n\tUSAGE: perl $0 <panam.ini file> \n\n\n\n";
# die "$USAGE" if ((scalar(@ARGV))< 1);
# my $option_file = $ARGV[0];
my $option_file = "panam.ini";
chomp($option_file);

die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";

my $usearch;
my $NGS_id_Results;
my $panam_output = "panam_output";
# my $parsing;
my $query_seq; 
my $user_file;
my $preprocess_output = "preprocess_output";
my $seq_F;
my $seq_R;
my $primF;
my $primR;
my $trim =0;
####
my $lib = "lib";
# my $fast = "$lib/FastTree"; # Tung: no need
# my $align = "$lib/hmmalign"; # Tung: no need
my $Reference = "Reference" ;
#my $build = "$lib/hmmbuild"; # Tung: no need
####
my $path = "bin/uclust3.0.617"; # change to bin 
my $us_version = "3.0.617";
my $dom;
my %dom;
my $aff_norm;


while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}
	
# get input options line
        else {
		
		if($line=~m/454_RUN_IDENTIFIER/gi) { 
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; 
			}
	 		die "\n\tOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( $NGS_id_Results ne "");
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "\n\tCould not create Directory $NGS_id_Results $!\n"; }
		}

# 		if (-d $NGS_id_Results."/".$panam_output) { `rm -r $NGS_id_Results/$panam_output` } 
		if (!(-d $NGS_id_Results."/".$panam_output)) { mkdir $NGS_id_Results."/".$panam_output || die "\n\tCould not create Directory $NGS_id_Results/$panam_output $!\n"; }


		if($line=~m/QUERY_SEQUENCES/gi) {
			$line=~s/QUERY_SEQUENCES//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi; 
			if($line) {
				$query_seq = $line; 
			}
			chomp ($query_seq);
			# See if the query sequences file exist and contains something
			die "\n\n\tCannot find query sequence file: $query_seq. Check panam.ini.\n\n" unless (-e $query_seq);
			die "\n\n\tQuery sequence file, $query_seq, appears to be empty!\n\n" if (-z $query_seq);
			
			my $ligne1= `sed -n '1p' $query_seq`;
			my $ligne2= `sed -n '2p' $query_seq`;	
			if (($ligne1 !~ /^>/) or ($ligne2 !~ /^[atcguATCGU]/)) {
				die "\n\n\tYour file $query_seq does not seem to contain sequences in fasta format.\n\n";
			}		
		}

		if($line=~m/USER_FILE/gi) {
			$line=~s/USER_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$user_file = $line; 
			}
			die "\n\n\t User file value missed. Check panam.ini.\n\n" unless (($user_file ne ""));
			die "\n\n\t Invalid user file value!\n\n" unless (($user_file eq "own") or ($user_file eq "pooledSamples") or ($user_file eq "eachSample"));	
		}

#		if($line=~m/REFERENCE_BASE/gi) {
#			$line=~s/REFERENCE_BASE//gi;
#			$line=~s/\t//gi;
#			$line=~s/^\s+//gi;
#			$line=~s/\s+$//gi;
#			if($line) {
#				$Reference = $line; 
#			}
#		}

		if($line=~m/FORWARD_PRIMER_NAME/gi) {
			$line=~s/FORWARD_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primF = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_NAME/gi) {
			$line=~s/REVERSE_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primR = $line; 
			}
		}

		if($line=~m/FORWARD_PRIMER_SEQUENCE/gi) {
			$line=~s/FORWARD_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_F = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_SEQUENCE/gi) {
			$line=~s/REVERSE_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_R = $line;
			}
		}

		############### normalisation
		if($line=~m/NBR_SEQ_NORM/gi) {
			$line=~s/NBR_SEQ_NORM//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$aff_norm = $line; #print "aff --- $aff_norm\n";
			} 
		}
		
		#####################################

		if($line=~m/DOMAIN/gi) {
			$line=~s/DOMAIN//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) { 
				$dom = $line;			
				if (defined $dom ) { #push (@dom, $dom);
					$dom{$dom} = 1
				}
			}
			foreach my $d (keys %dom) {
				die "\n\n\t Domain name is not correct. Check panam.ini.\n\n" unless (($d eq "bacteria") or ($d eq "eukaryota") or ($d eq "archaea"));
			}	
		}
	}
}

die "\n\n\t Invalid number of sequences to pick for the normalization! Check panam.ini.\n\n" if ((defined $aff_norm ) and ($aff_norm !~ /^\d+$/));
if (!(defined $aff_norm)) {$aff_norm = '200' }


print "$Reference\n";

############################################# trimmed profiles path ######################################################
# # # 
if (!(defined $primF) and (defined $primR))  {
	die "\n\n\t Primer forward not defined. Check panam.ini.\n\n";
}

if (defined $primF) {
	die "\n\n\t Primer reverse not defined. Check panam.ini.\n\n" unless (defined $primR);
	die "\n\n\t Primers' sequences missed. Check panam.ini.\n\n" unless ((defined $seq_F) and (defined $seq_R));
	die "\n\n\t Primers' sequences seem to be in the wrong format. Check panam.ini.\n\n" unless (($seq_F =~ /[ATCGUatcgu]/) and ($seq_R =~ /[ATCGUatcgu]/));
	$trim = 1
}

my @tax=();
my $Profiles;

if ($trim == 1) {
	#$Profiles = "Reference/Profiles_".$primF."_".$primR;
	$Profiles = "Reference/Trimmed_Profiles";
	if ((-d $Profiles)) {
		print "Trimmed profiles exist in $Profiles/\n";

	}
}

 
if ($trim == 0) {
	$Profiles= "Reference/Profiles";
}


###############################################################
########### hash query sequences 
my %h; my $c;
open (F, $query_seq);
while (<F>){
	if ($_ =~ /^>(.*?)\n/) {
		$c=$1;
		$c =~ s/\s+$//; 
		
	}
	else {
		chomp ($c);
		$h{$c}.=$_ ;
	}
}
close F;

########### hash tax
my %taxeuk;
open (T, "Reference/Taxonomy") || die "can not open file";
while (<T>) {
	my @e = split (/\t/, $_);
	chomp($e[0]); chomp ($e[1]); chomp ($e[2]);
	$taxeuk{$e[0]}=$e[1]."\t".$e[2];
}

#################################################
########### hash reference sequences + uclust
`mkdir $NGS_id_Results/$panam_output"/Similarity_annotation"` ;
`$path --sort $query_seq --output $NGS_id_Results/$panam_output/Similarity_annotation/seq_sorted 2>&1` ; 

print "usearch ...\n";

#my $repp = "Reference/Reference_all/Reference_bases/";
my $repp = "Reference/Reference_bases/";
my $e; my %seq;
#####comparer aux bases seed

my $rep_seed = $repp."/seed_bases/";
my @listeseed = <$rep_seed*_seed> ;
my %seed_affil;
foreach (@listeseed) {
	my $seed_file = $_;

	if ($seed_file =~ /\/Reference_bases\/seed_bases\/(.*?)_seed/ ){
		my $seed = $1;
		print "seed --- $seed\n";
		
		if (($us_version eq "1.1.579q") or ($us_version eq "3.0.617")) {
			`$path --input $NGS_id_Results/$panam_output/Similarity_annotation/seq_sorted_seed --lib $seed_file --uc $NGS_id_Results/$panam_output/Similarity_annotation/seq_results_seed_"$seed".uc --id 0.10 --libonly --allhits --maxaccepts 1 --rev `;
		}

		else {
			`$path --query $NGS_id_Results/$panam_output/Similarity_annotation/seq_sorted_seed --db $seed_file --uc $NGS_id_Results/$panam_output/Similarity_annotation/seq_results_seed_"$seed".uc  --id 0.10 --maxaccepts 1 --rev --local 2>&1`
		}	

		open (KK, $NGS_id_Results."/".$panam_output."/Similarity_annotation/seq_results_seed_".$seed.".uc") || die "can not open file" ;
		while (<KK>) {
			if (!($_=~ /\#/)) {
				my @t = split (/\t/, $_);
				if ($t[0] eq "H") {
					chomp($t[8]); chomp($t[3]); 
					$seed_affil{$t[8]}{$seed} = $t[3];
				}
			}
		} 
		close KK;		
	}
}

my %prim_affil ;
foreach my $u (keys %seed_affil) {
	my $max = 0;
	foreach my $uu (keys %{$seed_affil{$u}}) {
		if ($seed_affil{$u}{$uu} > $max) { 
			$max = $seed_affil{$u}{$uu}	; 		
			$prim_affil{$u} = $uu ;
		}
	}
}			

#open (DISC, ">".$NGS_id_Results."/".$panam_output."/Similarity_annotation/discarded.fasta") || die "can not open file" ;
foreach my $k (keys %prim_affil) {
	if (!(exists $dom{$prim_affil{$k}})) {
		#print DISC ">$k\n$h{$k}";
		delete ($h{$k});
	}
}
#close DISC;

open (QUERY, ">".$NGS_id_Results."/".$panam_output."/Similarity_annotation/query_seq") || die "can not open file" ;
foreach my $mm (keys %h) {
	print QUERY ">$mm\n$h{$mm}"
}		
close QUERY;

`$path --sort $NGS_id_Results/$panam_output/Similarity_annotation/query_seq --output $NGS_id_Results/$panam_output/Similarity_annotation/seq_sorted 2>&1` ; 
###################################################"" fin comparaison aux seed bases

my @listeficc = <$repp*.fasta>;
my $num_uc=1;
foreach(@listeficc) { 
	my $ficTemp = $_; 

	foreach my $domain (keys %dom) {
	if ($ficTemp =~ /$domain/) {

		#################################################
			if (($us_version eq "1.1.579q") or ($us_version eq "3.0.617")) {
				`$path --input $NGS_id_Results/$panam_output/Similarity_annotation/seq_sorted --lib $ficTemp --uc $NGS_id_Results/$panam_output/Similarity_annotation/seq_results_"$num_uc".uc --id 0.10 --libonly --allhits --maxaccepts 5 --rev `;
			}

			else {
				`$path --query $NGS_id_Results/$panam_output/Similarity_annotation/seq_sorted --db $ficTemp --uc $NGS_id_Results/$panam_output/Similarity_annotation/seq_results_"$num_uc".uc  --id 0.10 --maxaccepts 5 --rev --local 2>&1`
			}
		##################################################


		`cat $NGS_id_Results/$panam_output/Similarity_annotation/seq_results_"$num_uc".uc >> $NGS_id_Results/$panam_output/Similarity_annotation/seq_results.uc` ;
# 		
		$num_uc++;
		open(Q, "$ficTemp");
		while (<Q>){
			if ($_ =~ /^>(.*?)\s+/) {
				$e=$1;
			}
			else {
				chomp ($_);
				$seq{$e}.=$_ ;
			}
		}
		close Q ;
	}}
}


####################### sorting sequences
 
if (!(-z $NGS_id_Results."/".$panam_output."/Similarity_annotation/seq_results.uc")) {
	print L "usearch done\n";
} else {
	die("\#\#ABORT\#\#UCLUST_OUTPUT_EMPTY"); 
}
# 
print "Sorting sequences ...\n"; 

#my @alltax = ();
open (IL, "<Reference/sortingtax");
#@alltax=<IL>;
@tax=<IL>;
close IL;

#my @mydom = keys %dom;
#foreach my $grp (@alltax) {
#	chomp($grp);
#	my @parts = split("_",$grp);
#	my $type = pop(@parts);
#	if (grep $_ eq $type, @mydom) {# if (in_array($type, @mydom) ...
#		push(@tax,$grp); 
#	}
#}

open (F, $NGS_id_Results."/".$panam_output."/Similarity_annotation/seq_results.uc") || die "can not open file $NGS_id_Results/$panam_output/Similarity_annotation/seq_results.uc" ;
open (B, ">".$NGS_id_Results."/".$panam_output."/Similarity_annotation/best_hit_uc");
my %hits; my %target_hits;
while (<F>) {
	if (!($_=~ /\#/)) {
		my @t = split (/\t/, $_);
		if ($t[0] eq "H") {
			chomp ($t[8]) ,chomp($t[9]); chomp ($t[3]) ;
			$t[8] =~ s/\s+$//; $t[9] =~ s/\s+$//; $t[3] =~ s/\s+$//;
			$target_hits{$t[8]}{$t[9]} = $t[3];

			if ($t[4] eq "-") {
				my $rev = reverse ($h{$t[8]});
				$rev =~ tr/ACGTUacgtu/TGCAAtgcaa/ ; ##### modif 24/5
				chomp ($rev);
				$h{$t[8]} =$rev;
			}
		}
	}
}
 
my %tax_hits; my %seq_hits;
foreach my $query (keys %target_hits) {
	my $nb_target=0;
	foreach my $target (sort {$target_hits{$query}{$b} <=> $target_hits{$query}{$a}} keys %{$target_hits{$query}}) {
		my $ligne_taxo = $taxeuk{$target};
		my @t = split (/\t/, $ligne_taxo);
		print B "$query\t$target\t$t[0]\t$t[1]\n";
		$tax_hits{$query}{$nb_target} = $t[1];
		$seq_hits{$query}{$nb_target} = $target;
		$nb_target++;
		if ($nb_target == 5) {last}
	}
}

close B;


`mkdir $NGS_id_Results/$panam_output"/Similarity_annotation/Sorted_Sequences"`;
foreach my $tax (@tax) {
	chomp ($tax);
	my %tax_bis;
	if (-e $Profiles."/fasta/".$tax.".fasta") { 
		open (T, $Profiles."/fasta/".$tax.".fasta") || die "can not open file"; 
		my $aa; my %seq_ref;
		while (<T>) {
			if ($_=~ />(.*?)\s+/) {
				$aa = $1;
				chomp ($aa);
				$seq_ref{$aa}=1;
			}
		}
		close T;
		
		foreach my $seq_c (keys %seq_hits) {
# 			if ($tax_hits{$seq_c}{'0'} =~ /$tax/) {
			if ($tax =~ /$tax_hits{$seq_c}{'0'}/) {
				if (!(-d $NGS_id_Results."/".$panam_output."/Similarity_annotation/Sorted_Sequences/files_".$tax)) {
					`mkdir $NGS_id_Results/$panam_output"/Similarity_annotation/Sorted_Sequences/files_"$tax`;
				}
		
				open (TAX, ">>".$NGS_id_Results."/".$panam_output."/Similarity_annotation/Sorted_Sequences/files_".$tax."/seq_".$tax);
				print TAX ">$seq_c\n$h{$seq_c}\n";
		
				my %aze;
				foreach my $e (keys %{$seq_hits{$seq_c}}) {
					$aze{$seq_hits{$seq_c}{$e}}= 1
					
				}
	
				foreach my $uu (keys %seq_ref) {
					delete ($aze{$uu})
				}
	
				foreach my $qs (keys %aze) {
					$tax_bis{$qs} = 1;
				}
			}
		}
	
		foreach my $ww (keys %tax_bis) {
			print TAX ">$ww\n$seq{$ww}\n";
		}
		close TAX;
	}
}
